#import libraries
library(readr)
library(dplyr)
library(tidyr)
library(tidyverse)
library(ordinal)
library(effectsize)
library(emmeans)
library(lme4)
library(mediation)


#load wide-format data
data_wide <- read_csv("path", show_col_types = FALSE)
data_wide <- data_wide %>% mutate(id = row_number())
colnames(data_wide)


#transform the wide-format data to the long format
data_long <- data_wide %>%
  pivot_longer(
    cols = c(insurance_high, insurance_low, risk_high, risk_low),
    names_to  = c(".value", "inequality"),  # puts 'insurance' & 'risk' into value columns
    names_sep = "_"
  ) %>%
  rename(insurance_option = insurance) %>%   # optional: match your desired column name
  mutate(
    id               = factor(id),
    presentation_order            = factor(presentation_order, levels = c('low_inequality_first', 'high_inequality_first')),
    inequality       = factor(inequality, levels = c("low", "high")),
    insurance_option = factor(insurance_option)  # 3-level categorical outcome
  )
##transform insurance option to be continuous for ANOVA
data_long$insurance_option_num <- as.numeric(data_long$insurance_option) -1 

#get the sample information
##gender
gender_count <- table(data_wide$Gender) #male = 23, female = 26, non-binary or other = 1
gender_count 

##age
mean(data_wide$Age) #M = 41.80
sd(data_wide$Age) #SD = 11.47


#RM ANOVA on insurance choice
model_insurance <- aov(
  insurance_option_num ~ presentation_order * inequality + Error(id/inequality),
  data = data_long, 
)
summary(model_insurance)
##get effect size
effectsize:: eta_squared(model_insurance, partial = TRUE)
#######################################################
#Error: id
#                   Df Sum Sq Mean Sq F value Pr(>F)
#presentation_order  1   0.04  0.0371   0.051  0.823
#Residuals          48  35.20  0.7334               
#######################################################

#################################################################
#Error: id:inequality
#                             Df Sum Sq Mean Sq F value  Pr(>F)   
#inequality                     1  3.240   3.240  10.125 0.00257 ** (eta2 = .17)
#presentation_order:inequality  1  2.399   2.399   7.497 0.00864 ** (eta2 = .14)
#Residuals                     48 15.361   0.320                   
##################################################################################

##get descriptive
###get raw means on insurance choice
aggregate(insurance_option_num ~ inequality, data_long, mean) #low = .58, high = .94
aggregate(insurance_option_num ~ inequality, data_long, sd) #low = .67, high = .79

aggregate(insurance_option_num ~ inequality * presentation_order, data_long, mean)
aggregate(insurance_option_num ~ inequality * presentation_order, data_long, sd)
#########################################################
#  inequality    presentation_order insurance_option_num
#1        low  low_inequality_first            0.7407407 (sd = .76)
#2       high  low_inequality_first            0.8148148 (sd = .79)
#3        low high_inequality_first            0.3913043 (sd = .50)
#4       high high_inequality_first            1.0869565 (sd = .79)
########################################################

###get marginal means & contrasts
emmeans(model_insurance, ~ inequality * presentation_order)
####################################################################
#inequality presentation_order    emmean    SE   df lower.CL upper.CL
#low        low_inequality_first   0.742 0.144 81.3    0.456    1.028
#high       low_inequality_first   0.816 0.144 81.3    0.531    1.102
#low        high_inequality_first  0.393 0.147 85.1    0.100    0.686
#high       high_inequality_first  1.089 0.147 85.1    0.796    1.381
####################################################################

pairs(emmeans(model_insurance, ~ inequality * presentation_order))
###########################################################################################
#contrast                                               estimate    SE  df t.ratio p.value
#low low_inequality_first - high low_inequality_first    -0.0741 0.154 48.0  -0.481  0.9629
#low high_inequality_first - high high_inequality_first  -0.6957 0.167 48.0  -4.170  0.0007
###########################################################################################


#RM ANOVA on perceived financial risk
model_risk <- aov(
  risk ~ presentation_order * inequality + Error(id/inequality),
  data = data_long, 
)
summary(model_risk)
##get effect size
effectsize:: eta_squared(model_risk, partial = TRUE)
########################################################################
#Error: id:inequality
#                              Df Sum Sq Mean Sq F value  Pr(>F)   
#inequality                     1  18.49  18.490   9.335 0.00366 ** (eta2 = .16)
#presentation_order:inequality  1  11.94  11.938   6.027 0.01776 * (eta2 = .11)
#Residuals                     48  95.07   1.981       
########################################################################

####################################################################
#Error: id
#                   Df Sum Sq Mean Sq F value Pr(>F)
#presentation_order  1   1.05   1.047   0.375  0.543
#Residuals          48 134.20   2.796 
#####################################################################

##get descriptive
###get raw means on perceived financial risk
aggregate(risk ~ inequality, data_long, mean) #low = .2.92, high = .3.78
aggregate(risk ~ inequality, data_long, sd) #low = 1.50, high = 1.65

aggregate(risk ~ inequality * presentation_order, data_long, mean)
aggregate(risk ~ inequality * presentation_order, data_long, sd)
#######################################################
#inequality    presentation_order     risk
#1        low  low_inequality_first 3.333333 (sd = 1.36)
#2       high  low_inequality_first 3.555556 (sd = 1.45)
#3        low high_inequality_first 2.434783 (sd = 1.53)
#4       high high_inequality_first 4.043478 (sd = 1.85)
########################################################

emmeans(model_risk, ~ inequality * presentation_order)
#######################################################################
#inequality presentation_order    emmean    SE   df lower.CL upper.CL
#low        low_inequality_first    3.34 0.305 92.0     2.74     3.95
#high       low_inequality_first    3.56 0.305 92.0     2.96     4.17
#low        high_inequality_first   2.44 0.315 94.4     1.82     3.07
#high       high_inequality_first   4.05 0.315 94.4     3.43     4.68
#######################################################################

pairs(emmeans(model_risk, ~ inequality * presentation_order))
###############################################################################################
#contrast                                               estimate    SE   df t.ratio p.value
#low low_inequality_first - high low_inequality_first     -0.222 0.383 48.0  -0.580  0.9376
#low high_inequality_first - high high_inequality_first   -1.609 0.415 48.0  -3.876  0.0018
###############################################################################################


# Casual mediation analysis
data_long <- data_long %>%
  mutate(id = factor(id), inequality = factor(inequality)) %>%
  mutate(X = ifelse(inequality == levels(inequality)[1], -0.5, +0.5))

m.mod <- lmer(risk ~ X + (1 | id), data = data_long)
y.mod <- lmer(insurance_option_num ~ risk + X + (1 | id), data = data_long)

set.seed(123)
med <- mediate(m.mod, y.mod, treat = "X", mediator = "risk", sims = 10000, boot = FALSE)
summary(med)
